Packages & setup

# packages
if (!require("librarian")){
  install.packages("librarian")
  library(librarian)
}
librarian::shelf(
  calcofi/calcofi4r, 
  # marmap,
  cmocean, dplyr, DT, glue, here, interp,
  mapview, plotly, sf, skimr, tidyr)
mapviewOptions(fgb = F)

source(here("../apps/libs/db.R"))
d_vars    <- calcofi4r::get_variables()
d_cruises <- calcofi4r::get_cruises()
datatable(d_cruises)
# 1. choose cruise
(cruiseid <- d_cruises$cruiseid[1])
## [1] "2020-01-05-C-33RL"
# get casts, filtering by cruise
casts <- tbl(con, "ctd_casts") %>% 
  filter(cruiseid == !!cruiseid) %>% 
  select(cast_count, sta_id, date, longitude, latitude) %>% 
  collect() %>% 
  separate(
    sta_id, into = c("sta_line", "sta_offshore"), 
    sep = " ", convert = T, remove = F) %>% 
  mutate(
    day = difftime(date, min(date), units="days") %>% 
      as.integer()) %>% 
  st_as_sf(
    coords = c("longitude", "latitude"),
    crs = 4326)

mapview(casts, zcol="sta_line")
# mapview(casts, zcol="sta_offshore")
# mapview(casts, zcol="day")

table(casts$sta_line)
## 
##   60 63.3 66.7   70 73.3 76.7   80 81.7 81.8 83.3 85.4 86.7 86.8 88.5   90 91.7 
##    5    6    6    6    6    8    8    1    1   11    1   12    1    1   14    1 
## 93.3 93.4 
##   15    1
# remove station lines with only one cast
casts <- casts %>% 
  group_by(sta_line) %>% 
  mutate(sta_line_n = n()) %>% 
  filter(sta_line_n > 1)

table(casts$sta_line)
## 
##   60 63.3 66.7   70 73.3 76.7   80 83.3 86.7   90 93.3 
##    5    6    6    6    6    8    8   11   12   14   15
# 2. choose station line
(sta_line <- casts$sta_line[1])
## [1] 60
# filter by station line
casts <- casts %>% 
  filter(
    sta_line == !!sta_line)

bottles <- dbGetQuery(
  con,
  glue(
    "SELECT cast_count, depthm, t_degc 
    FROM ctd_bottles
    WHERE cast_count IN ({paste(casts$cast_count, collapse=',')})"))
d <- casts %>%
  inner_join(
    bottles,
    by="cast_count") %>% 
  st_drop_geometry()

d <- d %>% 
  select(sta_offshore, sta_line, depthm, t_degc) %>% 
  arrange(sta_offshore, sta_line, depthm)
datatable(d)

plotly

fig <- plot_ly(
  type='isosurface',
  x = d$sta_offshore,
  y = d$sta_line,
  z = d$depthm,
  value = d$t_degc,
  isomin = min(d$t_degc),
  isomax = max(d$t_degc))
fig
library(plotly)

fig <- plot_ly(
  type='isosurface',
  x = c(0,0,0,0,1,1,1,1),
  y = c(1,0,1,0,1,0,1,0),
  z = c(1,1,0,0,1,1,0,0),
  value = c(1,2,3,4,5,6,7,8),
  isomin=2,
  isomax=6
  )

fig

metR

shelf(metR)

d <- d %>% 
    mutate(depthm = depthm * -1)
  
ggplot(d, aes(sta_offshore, depthm, z = t_degc)) +
  geom_contour_fill(na.fill = TRUE)

ggplot(d, aes(sta_offshore, depthm)) +
  geom_contour_fill(aes(z = t_degc), kriging = T)

ggplot(d, aes(sta_offshore, depthm)) +
  geom_contour_fill(aes(z = t_degc), kriging = T) + 
  geom_point(size = 0.2)

ggplot2::geom_contour_filled()

# rever to positive depth
d <- d %>% 
    mutate(depthm = depthm * -1)
# View(d)
skim(d)
Data summary
Name d
Number of rows 99
Number of columns 4
_______________________
Column type frequency:
numeric 3
________________________
Group variables sta_line

Variable type: numeric

skim_variable sta_line n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
sta_offshore 60 0 1 73.13 12.58 53.00 60.00 70.00 80.00 90.0 ▇▁▆▆▆
depthm 60 0 1 157.24 154.88 0.00 30.00 100.00 250.00 516.0 ▇▃▂▁▁
t_degc 60 0 1 9.73 2.55 4.76 7.94 9.67 12.29 13.2 ▃▃▅▂▇
d %>% 
  ggplot(aes(x = sta_offshore, y = depthm)) +
  geom_contour_filled(aes(z = t_degc), bins = 10)+
  geom_point(size = 0.2) +
  scale_y_reverse() + 
  coord_cartesian(expand = F) +
  labs(x = "Offshore [id]", y = "Depth [m]")+
  theme_bw() %+%
  theme(
    panel.background  = element_rect(fill = "grey90"),
    panel.grid.major  = element_line(linetype = 3, colour = "grey60"),
    axis.text         = element_text(colour = 1, size = 10),
    axis.title        = element_text(colour = 1, size = 12),
    legend.background = element_blank(),
    legend.key        = element_blank(),
    legend.position   = "right") # +

  # metR::scale_x_longitude(ticks = 0.005,position = "bottom")
# Irregular data
# Use a dataset from the interp package
data(franke, package = "interp")

origdata <- as.data.frame(interp::franke.data(1, 1, franke))

grid <- with(origdata, interp::interp(x, y, z))

griddf <- subset(
  data.frame(
    x = rep(grid$x, nrow(grid$z)),
    y = rep(grid$y, each = ncol(grid$z)),
    z = as.numeric(grid$z)),
  !is.na(z))

ggplot(griddf, aes(x, y, z = z)) +
  geom_contour_filled() +
  geom_point(data = origdata)

  • cmocean: beautiful colormaps for oceanography
table(d$sta_offshore)
## 
## 53 60 70 80 90 
## 10 21 22 23 23
g <- with(
  d, 
  interp::interp(
    sta_offshore, depthm, t_degc))

gd <- subset(
  data.frame(
    x = rep(g$x, nrow(g$z)),
    y = rep(g$y, each = ncol(g$z)),
    z = as.numeric(g$z)),
  !is.na(z))

nbins = 8
cmocean_pal = "thermal"

p <- ggplot(gd, aes(x, y, z = z)) +
  geom_contour_filled(bins=nbins) +
  coord_cartesian(expand = F) +
  scale_y_reverse() + 
  labs(
    x     = "Offshore [id]", 
    y     = "Depth [m]",
    fill = "Temperature [C]") +
  scale_fill_manual(values = cmocean(cmocean_pal)(nbins)) +
  theme_bw()
p

# add points
p +
  geom_point(
    data = d %>%
      rename(
        x = sta_offshore,
        y = depthm,
        z = t_degc),
    size = 0.2)